%Plot figure 3

clear;
load ../Spectral/full_spectral_analysis.mat All Params Pall;

for ct=1:length(Pall); 
    Vexp(ct)=Pall(ct).Vexp;
    Vavg(ct)=mean(Pall(ct).V); 
    nw(ct)  =Params.nw(Pall(ct).setting(1));
    p(ct)   =Params.p(Pall(ct).setting(2));
    q(ct)   =Params.q(Pall(ct).setting(3));
    dt(ct)  =Params.dt(Pall(ct).setting(4));
    diff(ct)=Params.diff(Pall(ct).setting(5));
    S(ct)=All.S(ct);
end;

%specific formulation to plot
figure(1); clf; hold on; 

for ct=71, 
clf; hold on;
pl=find(All.S>4);
pl,
col='b';
clear P;
for ct2=1:length(pl); 
    P(ct2,:)=exp(Pall(ct).P(pl(ct2),:));
end;
s=Pall(ct).s';
Pm=nanmean(P);
Pm=Pm/nanmean(Pm);
s=s(6:end-2);
Pm=Pm(6:end-2); 
F=nanmean(Pall(ct).V)/Pall(ct).Vexp; 

%confidence intervals
k=length(pl)*2*nw(ct)-1;
dof=2*k; %degrees of freedom

%5 percent level
alpha=[0.95 0.05]; %alpha level 
conf1=gaminv(alpha,k,1/k); %5-95 confidence level
ptile1=[ones(size(s))*conf1(1); ones(size(s))*conf1(2)];

%1 percent level
alpha=[0.99 0.01]; %alpha level 
conf2=gaminv(alpha,k,1/k); %1-99 confidence level
ptile2=[ones(size(s))*conf2(1); ones(size(s))*conf2(2)];

%1 percent with Bonferroni, n=3
alpha=[1-0.01/3 0.01/3]; %alpha level 
conf3=gaminv(alpha,k,1/k); %1-99 confidence level
ptile3=[ones(size(s))*conf3(1); ones(size(s))*conf3(2)];

h=fill([s fliplr(s)],[ptile3(1,:) ones(size(s))],[0.8 0.8 0.8]);
set(h,'edgecolor',[0.8 0.8 0.8]);
h=fill([s fliplr(s)],[ptile2(1,:) ones(size(s))],[1 0.9 0.9]);
set(h,'edgecolor',[1 0.9 0.9]);
h=fill([s fliplr(s)],[ptile1(1,:) ones(size(s))],[1 0.75 0.75]);
set(h,'edgecolor',[1 0.75 0.75]);

plot(s,Pm,'k','linewidth',1);
axis tight;
h4=axis;
plot(1/41*[1 1],h4(3:4),'k--');
text(1/41,h4(3)+.02,' 1/41','fontsize',12);
font(gca,12);
h=xlabel('frequency (1/ky)'); font(h,14);
h=ylabel('normalized power density');  font(h,14);
title(num2str(ct)),
end;